/*************************************************************************
 *                                                                       *
 * This file is part of the FreeMAPs distribution.  For full copyright   *
 * information, see COPYRIGHT and COPYING.LESSER.                        *
 *                                                                       *
 * Author:           sqbang                                              *
 * Email:            sqbang@ustc.edu                                     *
 * Date:             2013.12.10                                          *
 *                                                                       *
 *************************************************************************/

//declaration header
#include "COMMON/tbox/MathUtilities.h"

//C++ headers
#include <cmath>
#include <float.h>
#include <limits.h>
#include <stdlib.h>

namespace FreeMAPs
{

template<>
bool MathUtilities<bool>::s_zero = false;

template<>
bool MathUtilities<bool>::s_one = true;

template<>
char MathUtilities<char>::s_zero = 0;

template<>
char MathUtilities<char>::s_one = 1;

template<>
int MathUtilities<int>::s_zero = 0;

template<>
int MathUtilities<int>::s_one = 1;

template<>
float MathUtilities<float>::s_zero = 0.0;

template<>
float MathUtilities<float>::s_one = 1.0;

template<>
double MathUtilities<double>::s_zero = 0.0;

template<>
double MathUtilities<double>::s_one = 1.0;

template<>
Complex MathUtilities<Complex>::s_zero = Complex(0.0, 0.0);

template<>
Complex MathUtilities<Complex>::s_one = Complex(1.0, 0.0);



template<>
bool MathUtilities<float>::is_NaN(
    const float& value)
{
    const int i = std::isnan(value);

    return (i != 0) ? true : false;
}



template<>
bool
MathUtilities<double>::is_NaN(
    const double& value)
{
    const int i = std::isnan(value);

    return (i != 0) ? true : false;
}



template<>
bool
MathUtilities<Complex>::is_NaN(
    const Complex& value)
{
    int i_re;
    int i_im;

    i_re = std::isnan(real(value));
    i_im = std::isnan(imag(value));

    return ((i_re != 0) || (i_im != 0)) ? true : false;
}



template<>
Complex
MathUtilities<Complex>::get_signaling_NaN()
{
    return Complex
           (
               std::numeric_limits<double>::signaling_NaN(),
               std::numeric_limits<double>::signaling_NaN()
           );
}



template<>
bool
MathUtilities<float>::equal_eps(
    const float& a,
    const float& b)
{
    float absmax = MathUtilities<float>::max(
                       MathUtilities<float>::abs(a),
                       MathUtilities<float>::abs(b));
    float numerator = MathUtilities<float>::abs(a - b);
    float denomenator =
        MathUtilities<float>::max(absmax, MathUtilities<float>::get_epsilon());

    return numerator / denomenator < sqrt(MathUtilities<float>::get_epsilon());
}



template<>
bool
MathUtilities<double>::equal_eps(
    const double& a,
    const double& b)
{
    double absmax = MathUtilities<double>::max(
                        MathUtilities<double>::abs(a),
                        MathUtilities<double>::abs(b));
    double numerator = MathUtilities<double>::abs(a - b);
    double denomenator =
        MathUtilities<double>::max(absmax, MathUtilities<double>::get_epsilon());

    return numerator / denomenator < sqrt(MathUtilities<double>::get_epsilon());
}



template<>
bool
MathUtilities<Complex>::equal_eps(
    const Complex& a,
    const Complex& b)
{
    double a_re = real(a);
    double a_im = imag(a);
    double b_re = real(b);
    double b_im = imag(b);

    return MathUtilities<double>::equal_eps(a_re, b_re) &&
           MathUtilities<double>::equal_eps(a_im, b_im);
}



template<>
Complex
MathUtilities<Complex>::min(
    Complex a,
    Complex b)
{
    return norm(a) < norm(b) ? a : b;
}



template<>
Complex
MathUtilities<Complex>::max(
    Complex a,
    Complex b)
{
    return norm(a) > norm(b) ? a : b;
}



template<>
int
MathUtilities<int>::abs(
    int a)
{
    return a > 0 ? a : -a;
}



template<>
float
MathUtilities<float>::abs(
    float a)
{
    return a > 0.0 ? a : -a;
}



template<>
double
MathUtilities<double>::abs(
    double a)
{
    return a > 0.0 ? a : -a;
}



template<class TYPE>
TYPE
round_internal(
    TYPE x)
{
    /* algorithm used from Steven G. Kargl */
    double t;
    if (x >= 0.0)
    {
        t = ceil(x);
        if (t - x > 0.5)
            t -= 1.0;
        return static_cast<TYPE>(t);
    }
    else
    {
        t = ceil(-x);
        if (t + x > 0.5)
            t -= 1.0;
        return static_cast<TYPE>(-t);
    }
}



template float round_internal<float>(float x);
template double round_internal<double>(double x);



template<>
float
MathUtilities<float>::round(
    float x)
{
    return round_internal<float>(x);
}



template<>
double
MathUtilities<double>::round(
    double x)
{
    return round_internal<double>(x);
}

}
